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Accreting black holes (BHs) produce intense radiation and powerful relativis- 
tic jets, which are affected by the BH's spin magnitude and direction. While 
thin disks might align with the BH spin axis via the Bardeen-Petterson ef- 
fect, this does not apply to jet systems with thick disks. We used fully three- 
dimensional general relativistic magnetohydrodynamical simulations to study 
accreting BHs with various BH spin vectors and disk thicknesses with mag- 
netic flux reaching saturation. Our simulations reveal a "magneto-spin align- 
ment" mechanism that causes magnetized disks and jets to align with the BH 
spin near BHs and further away to reorient with the outer disk. This mecha- 
nism has implications for the evolution of BH mass and spin, BH feedback on 
host galaxies, and resolved BH images for SgrA* and M87. 

Astrophysical black holes (BHs) operate as engines that convert gravitational binding energy 
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of accreting plasmas into intense radiation (1) and release BH spin energy (2-4) into powerful 
relativistic jets (5, 6). Relativistic jets from accreting BHs are commonly observed to emerge 
from active galactic nuclei (AGN) or quasars, x-ray binaries as microquasars, and gamma-ray 
burst (GRB) events. GRB jets allow one to probe the earliest epochs of star formation, whereas 
radiation and jets from AGN play a direct dynamical role via feedback that suppresses star 
formation in their host galaxies (7). 

BHs are also intrinsically interesting because they act as laboratories for probing Einstein's 
general relativity theory and for testing theories about accreting BHs and jets. Astrophysical 
BHs are characterized primarily by their mass (M) and dimensionless spin angular momentum 
(J). BHs have been measured to have masses of tens to billions of solar masses, like M87's BH 
with M ~ 6 x 10 9 M o (8). Spins have also been measured and span the whole range of possible 
values, including near the maximal value of j ~ 1 in the BH x-ray binary GRS 1915+105 with 
M ~ 14M Q and in the AGN MCG-6-30-15 with M ~ 3 x 1O 6 M (9, 10). Structures on BH event 
horizon length scales have recently been resolved by Earth- sized radio telescope interferometry 
for SgrA* (11, 12) and M87 (13, 14). 

The BH spin/rotational axis generally points in a direction that is tilted by some angle 6>tiit 
relative to the rotational axis of the plasma disk and its magnetic field's direction at large dis- 
tances. While the BH's present angular momentum axis is set by the history of plasma accretion 
and mergers with other BHs, the gas being currently supplied (at mass accretion rate M) to the 
BH can have an arbitrarily different angular momentum axis. This tilt influences the intensity 
of the radiation via changes in the gravitational potential felt by the plasma and also influences 
what an observer at different viewing angles sees due to disk warping and jet bending. 

One mechanism known to possibly affect the orientation of a disk or jet is the "Bardeen- 
Petterson" (BP) effect (15-18), where Lense-Thirring (LT) forces induced by the BH frame- 
dragging cause a misaligned disk to precess and warp until a local viscosity aligns a very thin 
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disk out to some distance (estimated to be out to r ~ I0r g -I0 5 r g , r g is a gravitational radius (19), 
depending upon assumptions) from the BH. Whereas the viscosity has been thought to result 
from turbulence driven by the magneto-rotational instability (MRI) that amplifies weak small- 
scale (< H, the disk height) magnetic fields (20), magnetohydrodynamical (MHD) simulations 
of weakly magnetized disks have not yet seen any BP alignment effect (21). The BP effect and 
LT precession remain commonly invoked mechanisms to understand the way tilt affects how 
BH mass and spin evolve (22, 23), how merging BHs are affected by any nearby plasma (24), 
and how disks and their jets are oriented (25, 26). 

Large-scale electromagnetic (EM) fields might also affect the jet's and disk's orientation via 
external confinement forces (27, 28). Estimates based upon large-scale magnetic fields being 
weaker than turbulent disk fields suggested that EM forces are insufficient to align the disk with 
the BH (27) or the BH with the disk (29, 30). Simulations without disks have given ambiguous 
results for the EM jet direction. For a uniform vertical magnetic field and no disk, the jet is 
directed along the magnetic field direction rather than the BH's tilted spin axis (31), whereas 
isolated magnetic threads tend to align with the BH spin axis when there is no disk to restrict 
their motion (32). 

We have used fully three-dimensional (3D) general relativistic (GR) MHD simulations (33) 
of accreting BHs to show that near the BH both the disk and jet reorient and align with the 
BH's spin axis. Our simulations were designed so that the magnetic field built-up to a natural 
saturation strength with, roughly, the disk's thermal+ram+gravitational forces balancing the 
disk's and jet's magnetic forces such that the trapped large-scale magnetic field threading the 
BH and disk became strong compared to the disk's turbulent field. The saturated field strength 
has been demonstrated to be independent of the strength of the initial magnetic field when 
the surrounding medium has a sufficient supply of magnetic flux, weakly dependent upon BH 
spin, and proportionally dependent upon disk thickness (4, 28, 34). We considered various BH 
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spins (35), BH tilts, and disks with a quasi-steady state height H to radius R ratio of H/R ~ 0.6 
for thick disks and H/R ~ 0.3 for thinner disks. Numerical convergence of our results was 
determined based upon both convergence quality measures for how well the MRI and turbulence 
was resolved (Table S2) and based upon explicit convergence testing (33). 

Let us motivate these MHD simulations by estimating whether EM forces are expected to 
dominate LT forces on the rotating heavy disk. Imagine a toy model with a flat heavy disk 
tilted and pushed-up against the magnetized jet generated directly by the rotating BH. The EM 
torque per unit area by the jet on the disk is r EM ~ rB r B^/A ~ r 2 B 2 Q. F /4 for magnetic field B 
bending on scale r with field line rotation frequency f2 F ~ j/8 for j ~ 1 (28). Introducing a 
dimensionless magnetic flux of Y « Q.l(Anr 2 B r )/ ^Anr 2 Mc with B in Gaussian units (r g and c 
reintroduced for dimensional clarity) that is consistent with measurements in our previous works 
(28), then r EM ~ T 2 MQ, F /(Snr 2 ). Meanwhile, the LT torque per unit area is t L t ~ ^lt^ with 
LT precession rate Q LT ~ 2j/r 3 , disk angular momentum per unit area L ~ Zri^, disk surface 
density 2 ~ M/(2nrv r ), and far from the horizon an effective viscosity a eS ~ v r /((H/R) 2 v^,). 
This gives t L t ~ jMv,p/(7rr 3 v r ). The ratio of the EM to LT torques for j ~ 1 is then t em /t L t ~ 
T 2 rv r /(64r g v $), with r g reintroduced for dimensional clarity. Far beyond the horizon, 

— — ~ -r 2r -^(H/R) 2 . (i) 

t lt 64 r g 

Over the horizon and in the jet, Y ~ 10 for our thinner disk models and Y ~ 17 for our thick disk 
models (28). Also, for both thicknesses, (r/r g )a eS ~ 15 and roughly constant with radius while 
Vr/v,/, ~ 1 near the horizon (28). So, at all distances, the jet's EM forces lead to t em /t L t ~ 2 
for our thinner disk models and t em /t L t k, 5 for thick disk models. So, we expect EM forces to 
dominate LT forces for both our thinner and thick disk models (including for small spins (33)). 

EM alignment forces are effective when they are larger than forces associated with the 
newly accreted rotating plasma with torque per unit area of r acc ~ Mv^/(2nr). So, r EM /r acc ~ 
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Y 2 Q F /(4rt^), and when these torques are equal one obtains an implicit equation for a "magneto- 
spin alignment" radius of 

£V g 2 T 2 

r msa ~ . j (2) 

(with r g reintroduced), within which EM forces can torque the accreting dense material. For 
sufficiently small Y or j, no alignment can occur. We obtain r msa > 30r g for our thinner and 
thick disk models that are sub-Keplerian by a factor 0.5 to 0.1, respectively (28), although 
accurate estimates require more physics (33) or simulations. 

Our self-consistent fully 3D GRMHD simulations started with a disk around an untilted BH 
where the BH spin axis, disk rotational axis, and emergent jet's direction all pointed in the verti- 
cal (z) direction. As the simulation proceeded, the mass and magnetic flux readily advected from 
large distances onto the BH. The magnetic flux vs. radius saturated on the BH and within the 
disk near the BH once magnetic forces balanced the disk's thermal +ram+gravitational forces. 
Magnetic braking causes such disks to become even more sub-Keplerian than weakly mag- 
netized thick disks (28), which means the classical thin disk inner-most stable circular orbit 
position is even less applicable than for weakly magnetized thick disks. The simulations were 
evolved for a long time period so that the disk reached a quasi-stationary magnetically saturated 
state out to about r ~ 40r g (28, 33). 

Then, the BH spin axis was instantly tilted by an angle of t ut,o ( see Table 1 for tilts used 
for different spins and disk thicknesses). The tilted disk-jet system underwent a violent rear- 
rangement for the larger tilts. The frame-dragging forces caused the nearly split-monopole BH 
magnetosphere to align with the BH spin axis, as expected because the misaligned angular mo- 
mentum was radiated away as part of the electromagnetic outflow on Alfven time scales. Then, 
the magnetic torque caused the heavy disk to lose its misaligned component of angular momen- 
tum and so reorient with the BH's rotating magnetosphere. The timescale for alignment seems 
to be roughly the Alfven crossing time near the heavy disk and BH. This "magneto-spin align- 
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ment" occurs because the magnetic field built-up to a natural saturation strength on the horizon 
where Y ~ 10 (depending upon the thickness and tilt), which led to the BH magneto sphere's 
forces dominating the disk dynamics and LT forces. 

All simulations (including with zero tilt) were then further evolved in time until all the 
tilted simulations reached a new quasi- stationary state out to r ~ 40r g . This ensured that any 
differences at later times in the disk and jet between tilted and untilted simulations were due 
to the BH tilt. This also ensured that each of the tilted and untilted simulations reached their 
own quasi-steady state values for the magnetic flux near the BH, magnetic flux in the disk, mass 
accretion rate, etc. We then measured the evolved relative tilt between the BH spin axis and 
the disk & jet at r = Ar g and r = 30r g (Table 1). For all our models, the disk & jet aligned 
with the BH spin axis near the BH, whereas the disk axis & jet direction deviated at larger 
distances. Such deviations are expected because the jet interacted with circulation with stronger 
mass inflows at larger distances (33). Despite the tilts and jet deviations, the BH's efficiency 
(defined as the ratio of energy out to energy in) was roughly 100% for j > 0.9 (Table SI), where 
more tilt led to reduced efficiency due to more spatially and temporally irregular mass inflow. 

Our most extreme case of a tilted BH accretion disk and jet system is the j = 0.99 model 
with a full tilt of # tilt , = 1.5708 * 90° and disk thickness H/R ~ 0.3. Even in this ex- 
tremely tilted case, the evolved disk and jet near the BH aligned with the BH spin axis (Fig. 1, 
Fig. SI, Movie SI). The jet's magnetic field winded around the persistent relativistic jet, and 
the magnetic field was well-ordered even for this highly tilted case. The jet was not symmetric 
around the jet axis, and instead there was a broad wing (with opening half-angle of about 25° 
by r = 40r g ) and a narrow wing (with opening half-angle of about 5° by r = 40r g ). At large 
distances from the BH, the jet drilled its way through the disk material and gradually got pushed 
away from the disk (Fig. 2 and Movie S2). The magnetic field winded around with a pitch angle 
of about 45° near the BH and had a smaller pitch angle at larger distances. By r ~ 300r g , the 
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Table 1 : Tilted Black Hole Disk- Jet Systems 



ModelName 


BH j 


Disk HjR 


Initial Tilt 


DiskTilt 

r=4r g 


Jet Tilt 

r=4r g 


DiskTilt 

r=30r g 


Jet Tilt 

r=30r g 


A0.94BfN40 


0.9375 


0.6 

















A0.94BfN40T0.35 


0.9375 


0.6 


0.35 


0.0 


0.0 


0.2 


0.2 


A0.94BfN40T0.7 


0.9375 


0.6 


0.70 


0.0 


0.0 


0.4 


0.3 


A0.94BfN40T1.5708 


0.9375 


0.6 


1.5708 


0.0 


0.1 


0.5 


0.7 


A-0.9N100 


-0.9 


0.3 

















A-0.9N100T0.15 


-0.9 


0.3 


0.15 


0.1 


0.1 


0.1 


0.2 


A-0.9N100T0.3 


-0.9 


0.3 


0.30 


0.2 


0.2 


0.2 


0.2 


A-0.9N100T0.6 


-0.9 


0.3 


0.60 


0.2 


0.3 


0.4 


0.3 


A-U.yJN 1UU1 1.3 /(Jo 


-u.y 


A O 

0.5 


1 C7AO 

1.3 /Uo 


A O 

0.1 


A A 

OA 


a n 

u.y 


A O 

U.o 


A0.9N100 


0.9 


0.3 

















A0.9N100T0.15 


0.9 


0.3 


0.15 


0.0 


0.0 


0.1 


0.1 


A0.9N100T0.3 


0.9 


0.3 


0.30 


0.1 


0.1 


0.2 


0.2 


A0.9N100T0.6 


0.9 


0.3 


0.60 


0.1 


0.1 


0.3 


0.3 


A0.9N10OT1.5708 


0.9 


0.3 


1.5708 


0.2 


0.3 


0.7 


0.6 


A0.9N100 


0.99 


0.3 

















A0.99N100T0.15 


0.99 


0.3 


0.15 


0.0 


0.1 


0.1 


0.1 


A0.99N100T0.3 


0.99 


0.3 


0.30 


0.1 


0.1 


0.2 


0.2 


A0.99N100T0.6 


0.99 


0.3 


0.60 


0.1 


0.1 


0.3 


0.4 


A0.99N100T1.5708 


0.99 


0.3 


1.5708 


0.1 


0.1 


0.6 


0.6 



Table 1: The simulation models listed by model name, which identifies the approximate BH 
spin of j = x (giving Ax), something about the magnetic field choices for some y and z (giving 
ByNz) (28), and the initial relative tilt (# t iit,o = P hi radians) between the BH spin axis and 
the disk rotation axis (giving Tp, where no T means 9^$ = 0). The second column gives 
the dimensionless BH spin (j), the third column gives the evolved quasi-steady state value of 
the disk height-to-radius ratio H/R between r ~ 20r g and r ~ 30r g , while the fourth column 
identifies the initial relative tilt (# t iit,o) between the disk+jet (having the same tilt initially) and 
the BH spin axis. The fifth and sixth columns give the evolved relative tilt between the BH spin 
axis and the disk & jet, respectively, at r = Ar g . The seventh and eighth columns give the same 
measurements at r = 30r g . A relative tilt of 0.0 means the disk or jet remained aligned with the 
BH spin axis, while a tilt equal to the initial tilt means the disk or jet was unaffected. 
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jet had become parallel with (but offset from) the outer disk rotational axis, and so the jet and 
counter-jet were also offset. 

Thus, our simulations have revealed a "magneto-spin alignment" mechanism that aligns the 
disk & jet axes with the BH spin axis near the BH once the magnetic field has saturated on the 
BH and within the disk (36). Unlike the BP effect, the mechanism actually works best for thick 
disks, and so the magneto-spin alignment mechanism should control jet systems where thick 
disks (due to accretion at either very low (37) or very high rates when H/R > 0.5) are expected. 
For example, for SgrA* and M87, if the BHs rotate sufficiently rapidly (33), then we expect 
their jet's and disk's photon spectra, temporal behaviors, and resolved images to be affected by 
non-zero relative tilts due to disk warping and jet bending near the BH. 

Tidal disruption flare events like Swift J 164449.3 +573451 are thought to be produced by 
very high accretion rates onto BHs, which launch fairly persistent jets that dissipate and give the 
observed emission (38). Our results suggest the inner disk and inner jet are both aligned with the 
BH spin axis, but the observed jet dissipating at large distances need not point along the BH spin 
axis. EM forces do not directly cause any precession (27), so the lack of LT precession-induced 
variability does not alone imply that the jet is necessarily driven by the BH spin power (26). 
Further, quasi-periodic oscillations (39) and long-term dips seen in this system's light curve 
might be explained by oscillations in the disk-jet magneto spheric interface (28) or by periods of 
magnetic flux accumulation and rejection by the BH (4) (both occurring for untilted systems) 
rather than by LT precession. 

Jet dissipation/emission (e.g. in blazars) might be due to the jet ramming into the disk until 
the jet aligns with the disk rotational axis at large distances. Measurements of BH spin in AGN 
and x-ray binaries might be affected by assumptions about the alignment between the disk, BH 
spin, and jet (9, 10). The cosmological evolution of BH mass and spin and AGN feedback for 
accretion at high rates might be affected by the higher BH spin-down rates and jet efficiencies 
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Figure 1: 3D snapshot for an evolved model with j = 0.99, initial relative tilt 6 tilt « 90°, and 
disk thickness of H/R ~ 0.3. The rotating BH sits at the center of the box of size r = -40r g 
to r = +40r g in each dimension. The snapshot shows the disk near the BH (yellow isosurface, 
which is mostly flat in the figure plane), the highly magnetized jet region (blue isosurface, 
with magnetic energy per unit rest-mass energy equal to about 70), the rotational axis of the 
disk both initially and at large distances (orange cylinder), outer disk (green-yellow volume 
rendering, more aligned with disk rotational axis at large distances), magnetic field vectors (like 
iron filings on that surface) for a cross-section of the jet (cyan vectors), and jet magnetic field 
lines (white lines) that trace from the BH out to large distances. The disk and jet near the BH 
are aligned with the BH spin axis and point mostly in-out of the figure plane, while at larger 
distances the jet points roughly half-way between the BH spin axis and the disk's rotational axis 
at large distances (pointing along the orange cylinder). 
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Figure 2: 3D snapshot that is similar to Fig. 1 but shows the outer disk (density isosurface 
in purple) at large distances from the BH in a box of size r = -350r g to r = 350r g in each 
dimension. The jet (blue isosurface, here with magnetic energy per unit rest-mass energy equal 
to about 4, still corresponding to the jet spine) is aligned with the BH spin near the BH but 
gradually gets pushed by the disk material and becomes parallel to (but offset from) the disk 
rotational axis at large distances. The strong interaction between the jet and disk has left an 
asymmetry or warp in the disk density at large radii. 
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compared to standard thin disk spin-down rates and radiative efficiencies (28) and also by how 
the jet aligns the disk material before LT torques can be effective so possibly leading to less 
change in the BH spin direction compared to the BP effect. 
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1 Methods 



In this supporting text, we describe our numerical methods, elaborate on details for the toy 
model in the report that estimates the effects of magnetic fields vs. Lense-Thirring precession 
around tilted black holes (BHs), present the equations of motion we numerically solve, discuss 
the diagnostics computed, and summarize the setup of our physical and numerical models. 

1.1 Numerical Methods 

We used fully three-dimensional (3D) general relativistic magnetohydrodynamic (GRMHD) 
simulations to study BH accretion flows with various BH spins (dimensionless value of J), 
relative tilts between the BH spin axis and the angular momentum (rotational) axis at large radii 
(#tiit,o m radians), and both thick and thinner disks (height H to radius R ratio of H/R ~ 0.6 and 
H/R ~ 0.3, respectively). 

Our simulations used the GRMHD code HARM based upon a conservative shock-capturing 
Godunov scheme with 2nd order Runge-Kutta time- stepping, Courant factor 0.8, Lax Friedrichs 
(LAXF) fluxes, simplified wave speeds, piecewise parabolic method (PPM) type interpolation 
(with no contact steepener, but with shock flattener based upon rest- mass flux density and total 
pressure), a staggered magnetic field representation, and a regular grid warping (40-42) as used 
in our prior works(28). A free version of HARM is available online for download (43). The 
HARM code has several key features that make it suitable for highly magnetized regions around 
BHs with tilted orientations (see appendix of (28)). In particular, the spherical polar coordinate 
axis is treated using transmissive (not reflective) boundary conditions so that the dissipation 
near the axis is kept to a minimum. This transmissive boundary condition treats all quantities 
in boundary cells near the axis as smooth and continuous functions whose values coincide with 
the non-boundary cells on the other side of the axis. This allows the flow and any components 
of the magnetic field to pass directly through the polar axis with little artificial dissipation while 
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maintaining a robust solution via the Godunov scheme. 

1.2 Toy Model 

In this section, we give more details about the toy model introduced in the report. 

The nature of jets around BHs has been a topic of interest for a long time (44-50). Our focus 
is on how the disk and jet are affected by a relative tilt with respect to the BH spin axis. Tilted 
BHs have been studied at an idealized level for a few decades. Solutions have been found for 
static scalar and vacuum magnetic fields around rotating BHs with relative tilt (51-55). For a 
vacuum magnetic field, it was found that the BH would be penetrated by non-zero magnetic flux 
for a maximally spinning (j = 1) BH only if there were relative tilt (54). Torque on a charged 
BH was also considered (56). For a review of these idealized works, see (55). Simulations have 
been required in order to push beyond these early idealized models, with some advances for 
simulations without tilt (for a review, see (57, 58)). Some axisymmetric simulations without tilt 
also see flow reversal near the BH (59), with the effect being quite strong once the magnetic field 
strength saturates (28). Untilted simulations have so far only studied how strong magnetic fields 
affect the aligned component of the BH spin (60) and (28). So far, no GRMHD simulations with 
tilt show any Bardeen-Petterson type alignment effect (21). Modern simulations of relatively 
thin weakly magnetized disks tilted relative to the BH spin axis have been studied, including the 
effect of tilt on the inner-most stable circular orbit (ISCO) (61), the development of shocks in 
the disk near the BH (62), the effect of tilt on observations of BH accretion disks (63, 64), and 
the effect on neutron star-BH mergers (65). Simulations of tilted disks with naturally saturated 
magnetic field strengths are required to understand how disks and jets are oriented under generic 
circumstances. Astrophysical applications justifiably only so-far invoke the Bardeen-Petterson 
effect as a possible alignment effect (48, 66-81), because no other type of alignment effect has 
yet been demonstrated in any GRMHD simulations and all theoretical arguments suggested that 



3 



magnetic torque alignment would not occur (27). 

Here we estimate the ratio of the electromagnetic (EM) force on the heavy disk relative to 
the Lense-Thirring (LT) force on the rotating heavy disk. This argument was used to get roughly 
correct the factors of order unity that appear in the report. We also account for some other subtle 
effects, like disk thinning, not accounted for in the report. Imagine a toy picture with a flat heavy 
disk tilted with respect to a rotating BH, with both the BH and disk threaded by a magnetic field 
(B) that bends on scale of radius r. The magnetic field is nearly co-rotating with the BH-disk 
system, which drives a current (J) into the ambient medium. Then, the EM torque per unit area 
on the disk by magnetic stresses is t E m = LEwi/dt ~ | J d9 sin Or x(JxB)| ~ 7rsmG^trB r B^/(4n) 
for radial field B r , toroidal field B^ ~ r sin 9Q, F B r for field line angular frequency of rotation of 
H F , and absolute magnetic flux on a radial shell of <D ~ 4nr 2 \B r \ with dimensionless magnetic 



flux of Y x 0.70/ yj4nrjMc for O in Gaussian units (28). We focus on the jet with a roughly 
constant magnetic flux O, while the magnetic field threading the disk would enhance the EM 
forces. The LT torque per unit area is r LT = dL LT /dt ~ |H LT x L| ~ sin # tilt Q LT L with |fi LT | = 
2|JbhI/t" 3 = 2j/r 3 for disk angular momentum per unit area L ~ ~Lrv^, disk surface density Z ~ 
M/(2nrv r ), effective viscosity a eS ~ v r /(G(H/R) 2 v^), and G ~ 1 accounts for GR corrections 
(28). Even untilted disks that reach magnetic flux saturation undergo magnetic braking that (in 
addition to the disk being thick) leads to a e ff ~ 1 and quite sub-Keplerian rotation (28), which 
makes us expect a) smaller LT forces ; b) EM forces can more easily change the disk rotational 
axis ; and c) the classical ISCO position is even less applicable than for weakly magnetized 
thick disks. The ratio of the EM to LT torques is then t em /t L t ~ T 2 (r g Q. ¥ )/(c j) (rv r )/(Sr g v^), 
with c and r g reintroduced for clarity. Far from the horizon, one can write this ratio as 



Assume the free parameters are j ~ 1, M{r) constant as valid for r < 30r g , Q, F ~ Q ff /4 = 





(1) 
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7 - /(8(l + V(l ~j 2 ))) ~ j/8 for the jet magneto sphere near the disk, and Y ~ 10-17 for our thinner 
and thick models, respectively (28). Also, over r ~ I2r g -20r g , the radially averaged a e s ~ 1- 
0.4 (with (r/r g )a e tf roughly constant over the solution) and H/R ~ 0.3-0.4, respectively for 
the thinner and thick models, where the thicker disk undergoes some compression by magnetic 
forces (28). For our thinner models we get t em /t L t ~ 2 while our thicker models get t em /t L t ~ 
4. So, we expect EM forces to dominate LT forces for roughly H/R > 0.2 at r ~ I0r g -20r g . 
In addition, near the horizon, v r /v$ ~ 1 (28), so near the horizon one obtains t em /t lt ~ 2 for 
thinner models and t em /t L t ~ 5 for thicker models. For non-zero spin, the value of T is roughly 
constant as a function of spin for both the thinner and thick disk models (28), so these estimates 
suggest EM torques dominate LT torques for all spins. 

The additional EM force by frame-dragging on the magnetic field threading the disk gives 
Y oc r in the disk (34)and naively Q F oc r^ 3 , which apparently leads to LT forces eventually 
dominating over the disk's internal EM forces at large radii (but the jet's EM forces external 
to the disk can still dominate LT forces as discussed above). However, significant BH spin 
energy and angular momentum is fed directly into the disk itself (28), which leads to angular 
momentum transport that forces the disk's magnetic field to rotate a bit closer to the BH mag- 
netosphere's value of Q, F even far from the BH (e.g. the rotation of the fluid flow completely 
re-aligns with the BH spin direction out to r ~ 40r g for untilted retrograde (j = -0.9375) thick 
disks (H/R ~ 0.6) such that the flow is no longer retrograde and has the same radial profile of 
but simply with an opposite sign compared to the otherwise identical model with j = +0.9375 
(28)). 

Using a similar analysis as above we can roughly estimate whether the jet's EM forces 
actually affect the flow. Such an analysis is easier than for the BP effect that requires a model of 
the assumed local viscosity that would be due to MHD turbulence, while large-scale EM fields 
directly affect the plasma. The torque per unit area associated with the newly accreted rotating 
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plasma at some fixed Eulerian radius is r acc ~ Mv^Klnr) so that r EM /T acc ~ T 2 Q. ¥ l(Arv ( j ) ). This 
gives a "magneto-spin alignment" radius of 



4vd 



(2) 



(with r g reintroduced) where the torques are equal, within which EM forces can torque the 
accreting dense material. For our thinner disk models and j ~ 1 with v$ oc r 1/2 (but sub- 
Keplerian in magnitude by a factor of 0.5) (28), the jet's external EM forces control the disk 
orientation out to r ~ 40r g . For our thick disk models and j ~ 1 with v$ oc r~ l (or Keplerian 
with a reduction factor of 0.1) out to r ~ 30r g up to where a steady-state was reached (28), EM 
forces can easily torque the disk out to r ~ 30r g and may even dominate out to r ~ 10000r g 
if the flow remains the same factor below Keplerian at such large radii (as is possible for thick 
disks that reach magnetic field strength saturation out to large radii). 

The magneto-spin alignment radius might depend upon spin. It is possible that for suffi- 
ciently slowly spinning BHs that produce weak jets, there might be no magneto-spin alignment 
effect. While rapidly rotating BHs might exist in SgrA* and M87 as some models/estimates 
show (82) and (14), this is still in debate (83, 84), so it is important to roughly estimate how spin 
might affect our results. For example, from Equation (2) and for the near-horizon region (i.e. 
r ~ r g ), horizon-scale alignment due to the jet alone occurs for j > 0.5 for our thinner disk mod- 
els and j > 0.2 for our thick disk models if we use ~ 0.5c as valid for the near-horizon when 
j ~ 1 (e.g., see figure 7 in (28)). However, more generally, over a range of spins, ~ r g Q, P near 
the horizon due to the BH dragging the magnetic field that then drags the plasma (4,34,28), so 
our estimates imply that all spins would allow for magneto-spin alignment of the disk close to 
the BH event horizon. 

The magneto-spin alignment radius might also depend upon the existence of mass inflow- 
outflow circulation that is commonly seen in thick disk simulations. The BH jet value of Y is 
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determined by force balance with the ingoing accretion near the horizon, but at larger radii one 
expects thick disks to produce winds that scale as M w i„d K r and so have an ingoing accretion 
rate following M in oc r (28), such that the EM forces are eventually ineffective against the heavy 
mass inflow forces. Stated another way, the local Y becomes too small in terms of the jet 
magnetic flux per unit square root of the local mass inflow rate. So, at some radius one expects 
the disk inflow to dominate the jet and for both the disk and jet to be parallel to the rotation axis 
of the outer disk (as consistent with the tilted simulations presented in this report). 

These estimates suggest that if an astrophysical system were to have a truncated thick disk 
inside r ~ r msa and still reach magnetic field strength saturation, then the entire disk and jet 
would mostly align with BH spin. 

These estimates also suggest that prior GRMHD simulations of tilted disks see no alignment 
effect (21)because their models have H/R ~ 0.1-0.2 and have a limited amount of magnetic 
flux available for accretion, both implying a small Y so the EM jet cannot dominate the disk 
dynamics and EM forces cannot dominate LT forces. For example, the standard "thick" torus 
with H/R ~ 0.2 with a single poloidal field loop in the initial torus is used in many GRMHD 
simulations. This setup is similar to the one we studied a few years ago (85) that we recently 
realized leads to limited magnetic flux with Y ~ 2 near the BH and so leads to low jet efficiencies 
despite a spin of j ~ 0.9 (28). Using Y ~ 2 in our above analysis shows that magneto-spin 
alignment does not occur for any normal range of spins (i.e. -1 < j < +1), which is consistent 
with existing tilted MHD simulations (21). 

Given the many non-trivial effects involved, simulations are optimal in order to obtain an 
accurate self-consistent solution once the magnetic field strength saturates so that all competing 
forces and effects can be accounted for. 
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1.3 Governing Equations 

We solved the GRMHD equations for a magnetized accretion flow around a tilted rotating BH 
defined by the Kerr metric in Kerr-Schild coordinates with internal coordinates x a = (t, jc (1) , jc (2) , jc (3) ) 
mapped to the spherical polar coordinates if = (t, r, 9, 4>). The z-axis of the spherical polar coor- 
dinate Kerr solution was rotated around the y-axis by #tiit,o (86). We write orthonormal vectors 
as u h contravariant (covariant) vectors as u l (w ; .), and higher-ranked coordinate basis tensors 
with no underbar. We work with Heaviside-Lorentz units, often set c = GM = 1, and let the 
horizon radius be r H . 

Mass conservation gives 

V^{p if) = 0, (3) 

where p is the rest-mass density, u M is the contravariant 4- velocity, and p = p u' is the lab-frame 
mass density. 

Energy-momentum conservation gives 

V? = o, (4) 

where the stress energy tensor Ty includes both matter (MA) and electromagnetic (EM) terms: 

T MA ^ = (po + e g +p g )uru v + p g ^ v , 
T EM ^ = b 2 ifu v + p h 8^-lfb v , 

The MA term can be decomposed into a particle (PA) term: T PAI ^ = pou v u M and an enthalpy 
(EN) term. The MA term can be reduced to a free thermo-kinetic energy (MAKE) term, which 
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is composed of free particle (PAKE) and enthalpy (EN) terms: 



= T^-p^Ja, (6) 

T EN t = (e g + Pg )y?u v + p g 6>l, 

such that r MAKE ^ = jpake^ + jen^ jj ere eg [ s fa e internal energy density and p g = (Y - \)e g 
is the ideal gas pressure with adiabatic index T = 4/3 (r = 5/3 may lead to somewhat different 
results ; (87, 88)). The sound speed is given as c s = ^Tp g /(p + e g + p g ). The contravari- 
ant fluid-frame magnetic 4-field is given by If, which is related to the lab-frame 3-field via 
If = E^hy/u 1 where hy = ifu v + 8y is a projection tensor that is the inverse metric on spa- 
tial hypersurfaces orthogonal to (giving an orthonormal comoving tetrad that is space-like 
diagonal) and that extracts the spatial part of the comoving magnetic field (giving b M u M = 0) 
and that in any basis gives W = (-m v 5 v )m /j + hyB y so the field has been split into parts parallel 
and perpendicular to it 1 , and 8y is the Kronecker delta function. The magnetic energy density 
(ub) and pressure (p h ) are u h = p h = lfb^/2 = b 2 12. The total pressure is p tot = p g + p b , 
and plasma /? p i asma = p g lpb- The 4- velocity of a zero angular momentum observer (ZAMO) is 
r] = {-a, 0, 0, 0} where a = 1/ y/-g" is the lapse, which is a convenient observer that follows 
worldlines with df/dt = -8' with B l = a 2 g tl such that rf = (l/a){l, -8'}. The BH metric can 
then be written as ds 2 = -a 2 dt 2 + h u (dx l + B l dt)(dx j + B j dt), which often simplifies calculations. 
We used this ZAMO to construct a frame that moves relative to the ZAMO with a time-like 
4-velocity of if = - yrf where y = -u a rj . Spatial interpolations of this relative 4-velocity 

— — — — — a 

always lead to time-like 4-velocities (unlike interpolations of the lab-frame 3-velocity) and also 
give a unique time-like observer even inside the BH ergosphere (unlike the lab-frame 4-velocity 
that can have two time-like observers inside the ergosphere). So this relative 4-velocity is well- 
suited to numerical simulations that rely heavily upon spatial interpolations and must not lose 



any information about the fluid frame in the ergosphere. For any 3-vector (e.g. B l ), the "quasi- 
orthonormal" vector is B t = B] yfgli computed in spherical polar coordinates. 
Magnetic flux conservation is given by the induction equation 

d t ( V=gff') = -dj[ V=g(5V - BV)], (V) 

where g = Det( ( g jUV ) is the metric's determinant, and the lab-frame 3-velocity is v l = u'/u'. No 
explicit viscosity or resistivity are included, but we used the energy conserving HARM scheme 
so all dissipation is captured (40, 89). 

The energy-momentum conservation equations are only modified due to so-called numerical 
density floors that keep the numerical code stable as described in detail in the Appendix of (28). 
The injected densities were tracked and removed from all calculations. 

1.4 Diagnostics 

The discussion in this section follows that of (28). Diagnostics were computed from snapshots 
produced every ~ 2r g /c. For quantities Q, averages over space «0) and time ([Q] t ) were 
performed directly on Q (e.g. on rather than on any intermediate values). Any flux ratio 
vs. time with numerator F N and denominator F D (F D often being mass or magnetic flux) was 
computed as R(t) = (F N (t))/[(F D )] t . Time-averages were then computed as [R] t . 

1.4.1 Fluxes and Averages 

For flux density F d , the flux integral is 

F(r) = J dA 23 F d , (8) 

where dA 2 3 = -y^gdx^dx^ (dA e(f) is the spherical polar version). For example, F d = pow (1) 
gives F = M, the rest-mass accretion rate. For weight w, the average of Q is 

f dA 0(P wQ 

Q w (r) = (Q) w = J - T -^ , (9) 

J dA e<p w 
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All 6, (p angles are integrated over. 



1.4.2 Disk Thickness Measurement 



The disk's geometric half-angular thickness is given by 



H/R,(((9-e f) p ) 



(10) 



where we integrated over all 6 for each r, <p, and 9 = n/2 + ((9 - n/2)) p was also integrated 
over all 9 for each r, 4>, and the final H/R(r) was from 0-averaging with no additional weight or 
Y 11 ^ factor. This way of forming H/R(r) works for slightly tilted thin disks or disordered thick 
disks. For a Gaussian distribution in density, this satisfies p/(p[9 = 0]) ~ exp(-9 2 /(2(H/R) 2 )). 

The value of H/R is computed at various locations and times, while for this work we report 
the value at the initial pressure maximum in the untilted solutions and at r ~ 30r g in the evolved 
quasi-steady state solutions. The thickness at large distances from the BH does not change 
significantly as the flow evolves (28)or after tilt is introduced in the simulation. Note that, unlike 
expected from modern thick disk accretion theory, the BH jet compresses the disk causing the 
geometric H/R to drop as the disk approaches the BH (28). 

1.4.3 Relative Tilt Measurement 

The tilt of the disk rotation axis and jet direction with respect to the BH spin axis was measured 
using the angular momentum (86). The BH angular momentum vector is 



with 0tii t! o is the initial tilt between the angular momentum vectors of the BH and the disk. The 
MHD angular momentum vector is 



7 B h = (aM sin 9 mfi x, 0, aM cos 9 m ,oz) , 



(ID 




(12) 



11 



in an asymptotically flat space-time, where 

(^mhd)p = 2 

for some frame with 4-velocity (for simplicity, we choose the Kerr-Schild lab-frame), and 

ZT = J dS a (xf T va - x v Tn , (14) 

where dL a = €p ySa dxPdy y dz 3 and the integral is performed for each radial shell. The unit vector 
y points along the axis around which the BH is tilted and the unit vector z points along the initial 
rotational axis of the disk. 

The plasma's relative tilt angle is(19) 

■^bh • AihdO") 



#tiit(r) = arccos 



(15) 



I-^bhIIAihdWI 

Just after the restart from a untilted simulation and instantly applying tilt, # t ii t (r) « tikO through- 
out the disk. 

As in (28), the disk is defined with the spatial integral performed with weight p u' (using p 
alone gives similar results) and normalized by the weight's volume integral. The jet is defined 
by taking the integral only over the region with b 2 /p > 1 with no weight. 

1.4.4 Efficiency of Energy Extraction Measurement 

The rest-mass flux, specific energy flux, and specific angular momentum flux are respectively 
given by 



M = ^Jp Q u r dA 6 



(16) 



E f TdA 9 # 

e=^- = -, (17) 

[M] t [M] t 

j J T r dA 6 t 

(18) 
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The net flow efficiency (efficiency of energy extraction from the BH) is given by 

E-M E EM (r) + E MAKE (r) 

Tj= : = : . (19) 

[Ml [M H ], 

Positive values correspond to an extraction of positive energy from the system at some radius. 
Values greater than 100% indicate that more energy goes out of the BH than goes in. 

1.4.5 Inflow Equilibrium and Quasi-Steady State 

Inflow equilibrium is defined as when the flow is in a complete quasi- steady- state and the accre- 
tion fluxes are constant (apart from noise) vs. radius and time. The inflow equilibrium timescale 
is 

'ie = N I dr I |, (20) 



I' 



for N inflow times from r = r ie and r, = 12r g to focus on the more self-similar part of the flow 
far from the BH. The value of r ie is measured directly from the simulation using Equation (20) 
rather than indirectly from an estimate of the a viscosity and disk scale-height. As in our prior 
work (28), we define the quasi-steady state out to some radius to mean that inflow equilibrium 
has been reached out to that radius under the condition that N is between 2 and 3. 

1.4.6 Modes and Correlation Lengths 

The flow structure was studied via the discrete Fourier transform of dq (related to quantity Q) 
along x = r,6,(f> giving amplitude a p for p = n,l, m, respectively. The averaged amplitude is 

N-l 



(\a p \) = (\T p (dq)\)= f 

J not x 



-2nipk 

, e^i" 

k=0 



Ydq> 



(21) 



computed at r = r H , 4r g , 8r g , 30r g . The x is one of r, 8, cf> and "not x" are others (e.g. 9, (p for 
x = r). The dq is (generally) a function of x on a uniform grid indexed by k of cells that 
span: 8r equal to 0.75r around r for x = r, n for x = 6, and 2n for x = <p. The Af is chosen so 
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all structure from the original grid is resolved, while the span covered allows many modes to be 
resolved. 

For all x,^ = ^f=gdx m dx {2) dx i3) 8Q/q N . For* = r, 9, we let q N = f nMx ^gdx (l) dx {2) dx {3) ([Q] t ), 
([Q]t) as the time-0 averaged Q, and 8Q = Q - ([Q] t )- Using dq removes gradients with r, 9 
so the Fourier transform acts on something closer to periodic with constant amplitude (see 
also (90)). For x = <p, we let q N = 1 and 6Q = Q because the equations of motion are (p- 
ignorable. For x = 9,(p, the radial integral is computed within +0.1 r. For x = r, 9, the (p integral 
is over all In. For x = r, <p, the 9 integral is over all n. For all x cases, the 9 range of values is 
over either the disk or jet (as defined earlier), where these regions are defined via 0-averaged 
quantities at each time. Notice we averaged the mode's absolute amplitude, because the am- 
plitude of (6Q) de -resolves power (e.g. m = 1 out of phase at different 9 gives (8Q) —> and 
a m —* 0) and is found to underestimate small-scale structure. 

We also computed the correlation length: A XtCor = x cor - x , where x = for x = 9, 4> and x 

is the inner radius of the above given radial span for x = r, where n COI = 6r/A rfior , l cor = n/A g ^ or , 

and m cor = (2n)/A^ cor . The Wiener- Khinchin theorem for the auto-correlation gives 

T~} x [<K>ol> 2 ] 
exp(-l) = , (22) 

^4[|<« P >ol> 2 ] 

where T~ l [{\a p> o\) 2 ] is the inverse discrete Fourier transform of (|a p |> 2 but with (a ) reset to 
(i.e. mean value is excluded). 

Turbulence (or flow structure in general) is resolved for grid cells per correlation length, 

Qp,cor = 7 ' (23) 

of Q P ,cor > 6, for x = r,9,(f> and p = n,l, m, respectively. Otherwise, modes are numerically 
damped on a dynamical timescale (even 2 = 5 would not indicate the mode is marginally 
resolved, because numerical noise can keep Q « 5 at increasing resolution until finally the 
mode is actually resolved - finally leading to an increasing Q > 6 with increasing resolution 
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; as seen by (91)). Reported Q PtC0I take 1/A X as the number of grid cells covering the span of 
A XtCor as centered on: middle of x m within the used radial span for x = r, 6 ~ n/2 (a reasonable 
approximation even for tilted models when this is measured at r « %r g ) for x = 9 for the disk 
and 6 ~ for x = 6 for the jet, and anywhere for x = <p. 

1.4.7 The Magneto-Rotational Instability 

The MRI is a linear instability with fastest growing wavelength of 

\v x a| 

I "rot I 

for x = 6,<p, where |v XjA | = sjb x b x /e is the ^-directed Alfven speed, e = b 2 + p + e g + p g , and 
rQ, rot = v mt . A M ri is accurate for Q. mt oc r~ 5/2 to r _1 . Q rot ,f A were separately angle-volume- 
averaged at each r, t. 

The MRI is resolved for grid cells per wavelength in the x = 9, <p directions, 

_ ^,mri 

of 2 > 10 and Q„ > 20 (92), where A r * dx m (dr/dx m ), A e * rdx {2) (d9/dx (2) ), and A * 
r sin 9dx^\d4>/dx i3) ). Volume-averaging was applied to t> x A /A x and |Q rot | were separately 6>, 0- 
volume- averaged before forming Q^mri- 

The MRI suppression factor corresponds to the number of MRI wavelengths across the full 
disk: 

c 2r(H/R) 

^MRI 

Wavelengths A < 0.5/i^MRi are stable, so the linear MRI is suppressed for S^mri < 1/2 when 
no unstable wavelengths fit within the full disk (93, 94). S^mri (or S^ wea k,MRi) uses averaging 
weight w = (b 2 p) 112 (or w = p), condition /? p i asma > 1, and excludes regions where density 
floors are activated. Weight w = (b 2 p) 1/2 is preferred, because much mass flows in current 
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sheets where the magnetic field vanishes and yet the MRI is irrelevant. When computing the 
averaged 5^ M ri, and |Q rot | were separately 6, <p- volume- averaged within +0.2r for each t, r. 
The averaged S^mri is at most 30% smaller than S</, W eak,MRi- 
A measure of the magnetic stress per unit pressure given by 

Pb 

can be used to determine how well-converged a numerical solution might be for MRI-dominated 
disk simulations. A value of a mag > 0.4 indicates a well-converged MRI-dominated simulation 
(92, 95). The numerator and denominator are separately volume averaged in 6, (f> for each r. Note 
that sin(20fc) = a mag for tilt angle 6 h (95). This measure is weakly applicable for our solutions 
where the MRI is suppressed and strong effects occur due to BH spin. However, at least for our 
thinner disk models this measure seems to roughly be applicable. For our thick disk models, 
the large amount of trapped magnetic flux and high BH spins lead to this measuring more about 
the angular momentum transported out from the BH into the disk than the MRI within the disk. 
For such thick disks, retrograde spins lead to negative a mag (28). 

1.5 Physical and Numerical Models 

This section describes our models. The model names are in the form AxByNzTp, where x is 
the approximate value of the BH spin, y identifies the field geometry (p=poloidal, f=flipping 
poloidal, t=toroidal), z identifies the normalization of the magnetic field, and p identifies the 
tilt (tilt is zero if not given). This is the same type of notation we have used previously(28), 
where for this paper only poloidal or poloidal flipping models are considered. For instance, our 
model A0.94BfN40T0.3 had a spinning BH (J = 0.9375), a poloidal field that flips polarity with 
radius, j8 p i asm a,min ~ 40 (i.e. smallest value of plasma /? p i asma is 40), and a tilt of 0.3 radians. 

All tilted models are continuations of the models without tilt given elsewhere(28)with pa- 
rameters given in their tables 1 and 2. The relative tilt was imposed (via a change in the BH spin 
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axis) at a time t = S000r g /c for thick disk models, t = lOOOOr^/c for thinner disk models except 
for the a = 0.99 models that used t = I5000r g /c. The thick disk models ran from 8000r g /c 
to roughly I7000r g /c, after which the initial magnetic field polarity inversion reaches the black 
hole, so no measurements are reported after this time. The thinner disk models with \ j\ = 0.9 
ran from lOOOOr^/c to roughly 22000r g /c, and the thinner disk models with j = 0.99 ran from 
15000r g /c to roughly 27000r g /c. These thinner disk models have no polarity inversions in the 
initial conditions. These durations allow the disk and jet to reach a quasi-steady state out to 
about r ~ 40r g within which results are reported. 

1.5.1 Physical Models 

The details of the original simulation's initial mass distribution (an isentropic hydro-equilibrium 
torus), magnetic field distribution (poloidal magnetic field loops), and velocity are as provided 
in (28). We only studied magnetic field configurations that led to a saturated magnetic field 
strength value near the BH, and we only studied cases that generate sufficiently large-scale 
magnetic field to produce a jet (85, 96). 

Our models were designed so that near the horizon the poloidal (R-z) magnetic flux reached 
a saturation point independent of the initial poloidal magnetic flux (4,34,28), which allowed us 
to provide robust results that are insensitive to any larger amount of magnetic flux available in 
the surrounding medium (97). In this sense, we identify the magnetic field strength reached as 
"natural" because it would require fine-tuning or an uncommon accident for the magnetic flux 
supply to be near saturation but somehow not simply far below or far beyond that required for 
near-BH saturation. A supply far below saturation would not lead to powerful jets as required 
by observations for many systems, while a large supply would simply lead to saturation out to 
some arbitrarily large distance and so not affect the near-BH dynamics. 

Our initial magnetic flux supply was still finite, and our H/R ~ 0.3 models did not (and 
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would not) reach a fully saturated state beyond r ~ 30r g (34), while our H/R ~ 0.6 models have 
enough magnetic flux supply to reach magnetic flux saturation out to r ~ 100r g (28). 

Our general relativistic magnetized disk simulations started with a tilted BH whose spin 
originally pointed towards the vertical (z) direction (28), and then (for this work) that BH spin 
direction was suddenly rotated by an angle 9^$ around another perpendicular axis (y-axis for 
definiteness) leading to the spin axis pointing more towards the jc-axis for larger # t nt,o- A value of 
#tiit,o = means no tilt, while 6^,0 = 90° means maximal tilt, while tilt;O = 180° means the BH 
was rotated so much that the spin value of j is simply the opposite sign. Hence, by symmetry 
we only need to consider # t iit,o = 0° to 90° for dimensionless spins in the range -1 < j < 1. 

For thick disks with height (H) to radius (R) ratio of H/R « 0.6 « 34° (the evolved quasi- 
steady state value of H/R at r ~ 30r g , where H/R « 0.6 at the initial torus pressure maximum 
at r = 100r g ), we obtained quasi-steady solutions with tilts # t iit,o = 0°,20°,40°, 90° with spin 
j = +0.9375. Note that j = -0.9375 is expected to behave similarly to a model with j = 
+0.9375 because these thick disk models result in the entire disk being forced to rotate in the 
same direction as the BH spin - giving a result that is completely independent of the initial disk 
angular momentum out to r ~ 40r g (28). 

For thinner disks with H/R ^ 0.3 « 17° (the evolved quasi-steady state value of H/R at 
r ~ 30r g , where H/R « 0.2 at the initial torus pressure maximum at r ~ 35r g ), we obtained 
quasi-steady solutions with tilts # tilti0 = 0°,9°, 17°,34°,90° with spins j = -0.9, +0.9, +0.99. 

Note that j is really a proxy for the actual BH angular frequency given by Q, H = jc/2r u 
(maximum value of 0.5) that sets the rate of frame-dragging. A value of j = 0.5, 0.9, 0.9375, 0.99, 1 
gives f2 H ~ 0.13,0.31,0.35,0.43,0.5. Hence, a value of j = 0.9 is really a middle-range value 
for the BH rotation rate. For small values of j < 0.5 that operate at less than a fifth of the 
maximum rotation rate, the effects of frame-dragging or tilt are expected to be none to minimal. 
Of course, in the limit of j = then the BH space-time is spherical and there is no physical 
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manifestation of any chosen # t iit,o- 

Table SI gives a list of diagnostics that we have measured as discussed in section 1.4. These 
diagnostics are measured just like in our prior work on untilted simulations (28). A list of iden- 
tical diagnostic measurements for other untilted simulations (including those used as a starting 
point for the tilted simulations) can be found in extensive tables elsewhere (28). The table shows 
that T and 77 remain relatively high for all tilted models, with some drop seen for fully tilted 
models. 

1.5.2 Numerical Models 

The resolution is chosen equal or higher than required to resolve turbulence and physical disk 
instabilities, as discussed previously (28). We chose a resolution N r x N g x that had a grid 
aspect ratio of 1:1:1 for most of the domain for r < 40r g . The grid is fully 3D with a 0-range 
of In. This allowed the <p dimension to be treated equally to the r - 6 dimensions. Thick disk 
models used a resolution of N r x N e x N<p = 272 x 128 x 256, while thinner disk models used 
N r x N e x N,p = 288 x 128 x 128 (twice larger Afy than our prior models for \j\ = 0.9). 

In our prior work on otherwise identical untilted simulations (28), we used numerous di- 
agnostics and numerical convergence quality measures to check that structures in the disk 
and jet were properly resolved (our measures are extensions of those discussed elsewhere 
(90-92, 95, 98)). In our prior study, for our thick disk models, we found that 272 x 128 x 256 was 
sufficient to overly resolve the MRI, well resolve the radial and azimuthal correlation lengths 
that describe the characteristic sizes of structures in the disk and jet, and marginally resolve the 
vertical correlation length due to magnetic Rayleigh-Taylor modes that lead to disk compres- 
sion. All thinner disk models were similarly resolved, except the \ j\ = 0.9 models that used 

= 64 did slightly under-resolve the azimuthal correlation length by a factor less than two. In 
convergence studies in that paper, the number of grid cells per azimuthal correlation length was 
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seen to scale with resolution as expected, so in this paper we doubled the resolution (from 
N$ = 64 to Nf = 128) to ensure that untilted thinner disk \ j\ = 0.9 simulations would properly 
resolve all correlation lengths. 

For these new tilted simulations, as in our prior work (28), we confirmed that the MRI and 
turbulence were resolved both initially and at late times through-out the disk and jet at many 
different radii. Table S2 gives a list of convergence quality measures as computed in section 1.4. 
These measures were computed the same way as in our prior work (28). A list of identical 
convergence quality measurements for the untilted simulations (used as a starting point for the 
tilted simulations) can be found in extensive tables elsewhere (28). Since the tilted simulations 
just continue the untilted simulations, the initial MRI convergence quality measures were the 
same as already reported elsewhere (28). In all the untilted models that we used for the tilted 
simulations, the initial MRI was well-resolved. 

In all cases, the convergence quality measures indicated that the simulations are at least 
marginally converged and some aspects were well-converged, which is expected given that was 
how we chose the resolutions as based upon similar measures for otherwise identical untilted 
simulations. 

As in our prior work (28), we also explicitly checked convergence of all our measures, 
including for disk tilt and jet tilt at late times, by repeating our thinner disk j = 0.99 t ut = 
1.5708 simulations with a resolution of 144 x 64 x 32 and 144 x 64 x 64. We did not double the 
resolution to 576 x 256 x 256 because such a single model would use about 16 million central 
processing unit (CPU)-hours on Kraken. Further, as we see next, convergence was already been 
obtained for T and tilt that we are focused on. 

As found in our prior work for untilted models, the value of T that controls the jet magnetic 
flux was insensitive to such changes in resolution. The disk and jet align with 6 rot « 0.1 at 
r « Ar g and 6 mt « 0.6 at r « 30r g as with the higher resolution model. The jet magnetic flux is 
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what leads to the EM force that aligns the disk with the BH, so the "magneto-spin alignment" 
effect is not sensitive to such changes in resolution. This is expected because the amount of 
trapped magnetic flux is set primarily by non-turbulent force balance against the large-scale 
magnetic field, and the EM jet power (and so force) does not strongly depend upon resolving 
the small-scale turbulence in the disk. 

The entire set of new tilted simulations (as continuations of untilted simulations) required 
about 22 million CPU-hours on the Extreme Science and Engineering Discovery Environment 
(XSEDE) Kraken cluster and another 11 million CPU-hours on the NAS Pleiades Westmere 
cluster. 
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Table SI: Simulation Physical Setup and Results 
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Table SI: The columns correspond to (1) the model name, (2) BH spin, (3) initial tilt, (4) final 
time of the simulation, (5) time- averaging range for any time-averaged quantity, (6) inner radius 
used to measure some quantities (<2 Xj mri,-> S^mru m Table S2), (7) outer radius used to measure 
some quantities (Q x ,mri , S dM Ri () in Table S2), (8) accurate disk thickness at r « 30r g , (9) mass 
accretion rate, (10) dimensionless magnetic flux on the horizon, and (11) BH efficiency. 

Movie SI 

This movie shows inner distances for the evolved state of the model with j = 0.99, t iit,o = 
1.5708, and H/R ~ 0.3 as shown in Figs. 1-2 in the report (See also http://youtu.be/XxIzazsE- 
sg ,20MB mp4 for full 1080p HD). The movie shows the full 3D structure that is difficult to 
see in a single 2D image because the BH, jet, and disk have axes that do not sit in the same 
two-dimensional plane. The movie shows the same rendering as Fig. 1 (without the vector field 
surface and with density volume rendering replaced by density isosurface) with the view rotated 
around a single axis so one can clearly see the full 3D structure of the BH spin direction, disk 
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Figure S 1 : Artistic 3D rendering showing an evolved 3D snapshot, similar to Fig. 1 in the report, 
for a model with j = 0.99, initial relative tilt # tilt « 90°, and disk thickness of H/R ~ 0.3. The 
rotating BH is shown as a black sphere and sits at the center of the box of size of about r = -40r g 
to r = +40r g in the vertical dimension and otherwise about r = -30r g to r = +30r g in the other 
dimensions. The jet is shown by the region where the magnetic energy density per unit rest- 
mass energy density is greater than about 10 up to about 100 (blue-white volume rendering). 
The logarithm of the density in the disk region is shown close to the BH (red-white volume 
rendering). The BH spin axis points vertically. Near the BH, the disk rotational axis and jet 
direction are strongly aligned with the BH spin axis. However, at t = and at large radius for 
all evolved times, the disk and jet axes pointed completely horizontally, i.e. perpendicular to 
the vertical axis. This shows the strong magnetic torques have caused the disk and jet to change 
orientation by 90° and to become aligned with the BH spin axis near the BH. At larger distances, 
the jet gradually bends away from the BH spin axis and aligns with the disk rotational axis at 
large distances. This leads to the jet and counter-jet being offset at large distances. 
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Table S2: Simulation Convergence Study 
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Table S2: The columns correspond to (1) the model name, (2) a mag with radial averaging over 
r = 12r g to r = 20r g , (3) number of grid cells per correlation length in the radial (n) direction, 
(4) number of grid cells per correlation length in the vertical-angular (/) direction, (5) number 
of grid cells per correlation length in the azimuthal (m) direction (all correlation lengths are 
measured for the rest-mass density and magnetic energy density at r « Sr g , where beyond the 
horizon all values are similar), (6) number of grid cells per vertical-angular MRI wavelength, 
(7) number of grid cells per azimuthal MRI wavelength, and (8) MRI suppression factor. The 
inner (?) and outer (o) radius values for columns 6,7,8 are measured at the inner and outer radii 
given in Table S 1 . 
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plane near the BH, jet direction, and disk rotational axis at larger distances (as represented by 
the orange cylinder). 

Movie S2 

This movie shows the outer distances for the evolved state of the model with j = 0.99, 9 iat ,o = 
1 .5708, and H/R ~ 0.3 as shown in Figs. 1-2 in the report (See also http://youtu.be/25WGtFVcHeO 
,11MB mp4 for full 1080p HD). The movie shows the same rendering as Fig. 2 with the view 
rotated around a single axis so one can clearly see the full 3D structure of the jet near the BH, 
the disk plane near the BH, the jet direction as it deviates at larger radii, the disk material that 
is pushed back into a warp by the jet, and the disk rotational axis at larger distances (as seen in 
the outer disk isosurface and represented by the orange cylinder). 
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